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ABSTRACT 

An aeroelastic analysis is developed which has general application to all 
types of axial-flow turbomachinery blades. The approach is based on lin- 
ear modal analysis, where the blade's dynamic response is represented 
as a linear combination of contributions from each of its in-vacuum free 
vibrational modes. A compressible linearized unsteady potential theory 
is used to model the flow over the oscillating blades. The two- 
dimensional unsteady flow is evaluated along several stacked 
axisymmetric strips along the span of the airfoil. The unsteady pressures 
at the blade surface are integrated to result in the generalized force acting 
on the blade due to simple harmonic motions. The unsteady aerodynamic 
forces are coupled to the blade normal modes in the frequency domain 
using modal analysis. An iterative eigenvalue problem is solved to de- 
termine the stability of the blade when the unsteady aerodynamic forces 
are included in the analysis. The approach is demonstrated by applying 
it to a high-energy subsonic turbine blade from a rocket engine 
turbopump power turbine. The results indicate that this turbine could 
undergo flutter in an edgewise mode of vibration. 


ACKNOWLEDGEMENTS 


I wish to thank my employer, Sverdrup Technology, Inc. for its continued 
support and encouragement during the course of my graduate study. In 
addition, I also would like to thank the NASA Lewis Research Center 
(LeRC), in particular George Stefko and Oral Mehmed, for their sponsor- 
ship of this research activity. 

On a personal note, I would like to thank Dr. Robert Kielb of General 
Electric Corp. for introducing me to the field of aeroelasticity. I would also 
like to express my warm regards to my friend Dr. Durbha Murthy, a resi- 
dent research associate at NASA LeRC, for the many hours of thoughtful 
discussions we shared over several years. 


ii 



TABLE OF CONTENTS 


ACKNOWLEDGEMENTS 

LIST OF FIGURES 

LIST OF TABLES 

NOMENCLATURE 

CHAPTER 1 - INTRODUCTION 

1.1 Introduction 

1 .2 Literature Review 

1.21 Turbomachinery Aeroelastic Problems 

1.22 Unsteady Aerodynamic Models .... 

1.23 Aeroelastic Analysis Models 

1.3 Objective of Current Work 

CHAPTER 2 - ANALYTICAL DEVELOPMENT . 
2.1 Aerodynamic Analysis 

2.11 Steady Full Potential Flow 

2.12 Linearized Unsteady Potential Flow . 

2.2 Structural Dynamic Analysis 


. ii 
vii 


ix 

x 


. 1 
. 1 

. 4 
. 5 
. 9 
16 
22 

24 

24 

25 
28 

33 


2.3 Aeroelastic Model . 


38 



CHAPTER 3 - APPLICATION OF METHOD 48 

3.1 Computer Program 48 

3.2 SSME HPOTP Turbine Description 49 

3.3 Aeroelastic Model 52 

3.4 Structural Dynamic Results 54 

3.5 Aerodynamic Results 56 

3.51 Steady Aerodynamic Results 58 

3.52 Unsteady Aerodynamic Results 59 

3.6 Aeroelastic Results 62 

CHAPTER 4 - CONCLUSION AND RECOMMENDATIONS 67 

4.1 Conclusion 67 

4.2 Recommendations 69 

REFERENCES 72 

APPENDIX A 77 

TABLES 88 

FIGURES 90 


IV 



LIST OF FIGURES 


Figure 1. 
Figure 2. 
Figure 3. 

Figure 4. 

Figure 5. 
Figure 6. 
Figure 7. 
Figure 8. 
Figure 9. 

Figure 10. 
Figure 1 1 . 
Figure 12. 
Figure 13. 
Figure 14. 
Figure 15. 
Figure 16. 
Figure 17. 


Cascade representation of rotor inlet and exit flows 
Cascade and airfoil geometry and nomenclature 
Steady solution blade-to-blade passage and boundary con- 
ditions 

Steady flow computational meshes, (a) global mesh, (b) local 
mesh 

Unsteady solution cascade boundary conditions 

Oscillating airfoil displacement vector r definition 

FREPS program aeroelastic stability logic flowchart 

Modal aerodynamic matrix [A] logic flowchart 

SSME high pressure oxygen turbopump (HPOTP) cross- 

section 

SSME HPOTP turbine blade finite element model 
Turbine blade aeroelastic strip definitions 
Turbine blade airfoil description at 75% span 
Calculated natural frequencies versus rotational speed 
Calculated deformed mode shapes for modes 1 and 2 
Calculated deformed mode shapes for modes 3 and 4 
Airfoil mode shape cross-sections at strip no. 1 
Turbine blade fluid velocity diagram at strip no. 1 


v 


Figure 18. 
Figure 19. 
Figure 20. 
Figure 21. 
Figure 22. 

Figure 23. 

Figure 24. 

Figure 25. 

Figure 26. 

Figure 27. 

Figure 28. 


Steady flow Mach number contours for strip no. 1 
Steady flow Mach number contours for strip no. 3 
Turbine blade surface Mach number distributions 
Turbine blade surface static pressure distributions 
Turbine blade unsteady surface pressure distributions (real 
part) due to motion in second mode (edgewise) 

Turbine blade unsteady surface pressure distributions (im- 
aginary part) due to motion in second mode (edgewise) 
Turbine blade unsteady aerodynamic work per cycle distrib- 
ution due to motion in second mode (edgewise) 

HPOTP turbine rotor tuned root locus for vibration in first 
mode (bending) 

HPOTP turbine rotor tuned root locus for vibration in third 

mode (edgewise), mechanical damping neglected 

HPOTP turbine rotor tuned root locus for vibration in second 

mode (edgewise), mechanical damping neglected 

HPOTP turbine rotor tuned root locus for vibration in second 

mode (edgewise), 1 % mechanical damping included 


VI 



LIST OF TABLES 


Table I. HPOTP turbine aeroelastic strip properties 
Table II. HPOTP turbine calculated natural frequencies 


NOMENCLATURE 


Symbols 

a 

[A] 

b 

S* 

c 

C p 

C P 

[C] 

rc G j = 


= Fluid sonic velocity 


= Modal unsteady aerodynamic matrix, eq (23) 


Airfoil semichord length 


c cos £ Airfoil projected chord length at stagger angle 


2b Airfoil chord length 



Static pressure coefficient 


P 

pV 2 


Harmonic unsteady pressure coefficient 


Mechanical and structural damping matrix 


Generalized damping matrix, eq (14) 


viii 


{F A (t)} = External aerodynamic forcing function 


{F M {t)} = External mechanical forcing function 


h = fluid specific enthalpy 




k = Reduced frequency 


[K] = Structural stiffness matrix 


£K g 3 = Generalized stiffness matrix, eq (14) 


[K‘] = state-space modified stiffness matrix, eq (27) 


M = Mach number 


[/W] = Structural mass matrix 


= Generalized mass matrix, eq (14) 


IX 


= State-space modified mass matrix, eq (27) 


[M*] 

n = Number of degrees-of-freedom in physical space 

n = Airfoil surface normal vector 

N = Number of rotor blades 

P = Fluid pressure 

P = Steady-state fluid pressure 

P = Harmonic unsteady fluid pressure 

{q} = Vector of time dependent modal coordinates 

{q} = Vector of harmonic modal coordinates, eq (19) 

{0} = Vector of generalized forces, eq (15) 

{Q*} = Vector of aerodynamic generalized forces, eq (16) 


x 



{Qm} = Vector of mechanical generalized forces, eq ( 16 ) 


r = Airfoil displacement vector 


{R} = Vector of modal unsteady aerodynamic forces 


[/?] = Modal unsteady aerodynamic force matrix 


t = time 


{u} = Vector of physical displacements 


V = Fluid velocity vector 


l/ T = Fluid velocity vector tangent to airfoil surface 


xi 


Greek 


a = Absolute flow angle 


P = Relative flow angle 


= Gas specific heat ratio 


= Modal damping ratio 


X = Complex aeroelastic eigenvalue, eq (30) 


M = Real part of complex eigenvalue X, eq (30) 


v = Imaginary part of complex eigenvalue X, eq (30) 


£ = Cascade stagger angle measured from axial direction 


{£} = Vector of time-dependent state-space modal coordinates, eq 

(25) 


{£} = Vector of harmonic state-space modal coordinates, eq (25) 


p 


= Fluid density 


o = Interblade phase angle 


= Airfoil surface tangent vector 


<f> = Scalar velocity potential 


<f> = Harmonic unsteady velocity potential 


ct> = Steady-state velocity potential 


{</>„} = Free-vibration natural mode shape vector 


[O] = Free-vibration modal matrix 


co = Radian frequency 


cot = Assumed flutter frequency 


co n = Free-vibration natural frequency 




CHAPTER 1 - INTRODUCTION 


1.1 Introduction 


The development of propulsion system technology during the last forty 
years has encountered and overcome many technological barriers. Se- 
veral problems associated with advanced high temperature materials, 
turbine cooling, and fuels and combustion were resolved with the end re- 
sult of significantly higher component efficiencies and reduced fuel con- 
sumption. For gas turbine powerplants, these advances led to lighter 
overall designs and higher power densities compared to earlier designs. 
The accomplishments of lighter designs for the turbomachinery compo- 
nents also Jed to some drawbacks due to the reduced margins on the de- 
sign facto r-of-safety. 

Fan and compressor designs during this period placed primary emphasis 
on higher aerodynamic loadings, and lighter blade weight to increase the 
gross power without sacrificing fuel efficiency. This trend led to increas- 
ingly frequent problems resulting from aerodynamic and structural inter- 
actions. A particularly serious problem was due to aeroelastic instability, 
flutter, which began to occur more often. This was a surprising problem 

because the technical consensus at the time considered that 
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turbomachinery blading had sufficiently high stiffness and strength to 
avoid aeroelastic problems. 

Aeroelasticity is a science concerned with the mutual interactions among 
the structural (inertial and elastic) and aerodynamic characteristics of a 
structure immersed in a flowing fluid. Under certain conditions, when a 
vibrating structure begins to extract energy from the flowing fluid the 
structure will experience self-excited vibration. This self-excited vibration 
is a dynamically unstable condition, and is referred to as flutter. Up to 
1960, the most common use of aeroelasticity theory was focused on the 
analysis of aircraft structures. Static aeroelasticity associated with wing 
divergence and dynamic aeroelasticity and flutter control were areas of 
active research for structures subjected to external flow, such as aircraft 
structures. 

The aeroelastic problems which began to occur for systems under internal 
flow, i.e. turbomachinery, resulted from a combination of (a) lighter and 
higher stressed blade designs operating at (b) higher fluid velocities rel- 
ative to fixed wing external aircraft structures. Successful solution of the 
aeroelastic problems in turbomachinery required more robust mechanical 
designs like the incorporation of part-span and tip shrouds and the use 
of lacing wires and intentionally introduced friction damping. Unfortu- 
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nately, these mechanical modifications resulted in an increase in compo- 
nent weight in addition to a degradation in the aerodynamic performance 
of the machine. 

Designers of modern propulsion systems have discovered that multi- 
disciplinary effects must be considered in order to achieve reliable 
turbomachinery designs. The interaction of the distinct physical proc- 
esses of fluid dynamics, structural dynamics and thermodynamics of the 
entire engine system is becoming more important. The traditional method 
of designing for aerodynamic performance, independently of the struc- 
tural dynamics of the system is no longer a viable approach. 

The application of aeroelasticity theory to propulsion systems has enjoyed 
considerable attention during the time period from 1960 to the present. 
Advances in methods for computational fluid dynamics (CFD) and finite 
element analysis methods have greatly expanded the applications of 
aeroelastic theory. Current methods for predicting aeroelastic stability 
(flutter) rely upon a variety of unsteady aerodynamic models, from simple 
flat-plate representations of the blade up to complicated CFD analyses. 
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1.2 Literature Review 

A thorough review of turbomachinery aeroelasticity requires an under- 
standing of the physical mechanisms which dominate this type of prob- 
lem. There are basically two physical processes which interact with one 
another within an aeroelastic system. These processes are (a) the un- 
steady aerodynamic behavior and (b) the structural dynamic behavior of 
the turbomachinery blading. The unsteady aerodynamic problem is con- 
cerned with estimating the aerodynamic response of a blade which is 
undergoing a periodic oscillation. This problem is complicated due to the 
fact that there is a transfer of energy between the blade and the flowing 
fluid. The structural dynamics of turbomachinery blading requires the 
study of the dynamic response of the blade under both free and forced 
vibration conditions. Aeroelasticity is an investigation of these combined 
aerodynamic and structural dynamic effects. 

The review of previous work within the field of turbomachinery 
aeroelasticity is presented in the following order. A survey of the occur- 
rence of aeroelastic problems, most notably self-excited vibrations or 
flutter, for various types of turbomachinery including aircraft engine fans, 
compressors and turbines is presented. This is followed by an examina- 
tion of some of the fluid dynamic models which have been proposed to 
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solve the unsteady aerodynamic problem of oscillating airfoils under se- 
veral different flow regimes. Finally, a review of aeroelastic models which 
have been applied to turbomachinery analyses over the past thirty years 
is included. 

1.21 Turbomachinery Aeroelastic Problems 

Reports of aeroelastic instability, or flutter, occurring in actual 
turbomachinery such as aircraft engines or power generation turbines are 
relatively rare. This is not unexpected, because the design and develop- 
ment of high-energy gas turbines represents a very expensive undertak- 
ing, often requiring millions of dollars and many years of design work 
prior to customer acceptance or certification. Some instabilities are re- 
ported in the literature, either as a result of catastrophic failures, or de- 
velopment delays in highly visible aerospace projects. Gas turbine 
manufacturers are often quick to publish results when such problems 
come into common knowledge in order to prove that a redesign will per- 
form successfully. 

Carter and Kilpatrick [1] were one of the first to publish an investigation 
of the aeroelastic response of a multi-stage compressor operating under 
part-speed conditions. Their experiments were conducted on the 
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compressors of power generation turbines to determine the vibratory 
stresses on the stator vanes during stall flutter and rotating stall condi- 
tions. The results from their experiments led to the introduction of some 
modified design correlations based on reduced frequency and incidence 
angle limits for high-speed compressors. 

Jeffers and Meece [2] present a diagnosis of a fan stall flutter occurrence 
for the Pratt & Whitney FI 00 turbofan engine which was discovered during 
engine development testing. Stall flutter was found to occur for the first- 
stage fan at high Mach numbers under off-design conditions near the 
surge threshold of the engine. The fan was designed using the standard 
empirical correlation rules employed for all of Pratt & Whitney's engine 
designs up until that year. A full structural redesign of the fan including 
airfoil thickening and the modification of the part-span shroud was re- 
quired to overcome the flutter problem. The authors expressed the need 
for a more reliable subsonic flutter prediction system to augment the ex- 
isting empirical database which at that time was the standard technique 
used for aeroelastic design. 

An outline of the aeroelastic design system used by General Electric was 
described in the work of Cardinale, Bankhead and McKay [3]. This paper 
presents a general overview of the empirical design system which had 
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been verified over several years of engine testing and in-service opera- 
tion of aircraft gas turbine engines. A detailed discussion of 
turbomachinery mechanical and aerodynamic design, including the ex- 
perimental identification of aeroelastic problems is also included. A 
strong conclusion from this paper is that the inherent complexities of 
turbomachinery, due to variations in operating conditions and flight en- 
velopes, require the application of conservative design rules based on 
thousands of hours of existing engine test experience. 

Aeroelastic problems are not solely limited to aircraft gas turbine engines. 
However, these engines are more prone to these problems because of the 
stringent thrust/weight requirements of propulsion systems relative to 
other types of gas turbines. A flutter problem encountered within the last 
stage turbine of a Westinghouse industrial combustion turbine was docu- 
mented by Scalzo, Allen, and Antos [4]. The problem was found to occur 
within the last stage turbine under operating conditions where the inlet 
air density was high relative to normal operation. High-cycle fatigue 
cracking was caused by a self-excited vibration in the first bending mode 
of the blade. A rotor-based telemetry system was used to quantify the 
vibratory stress levels under hostile high temperature operating condi- 
tions [5], The turbine blade was redesigned to stiffen the base of the 
blade and to increase the natural frequency of the blade, thus giving more 
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acceptable margin for flutter. Large amplitude dynamic response due to 
buffeting was also reported for the blades which did not flutter, although 
the large response is probably attributable to extremely low aerodynamic 
damping for those blades. 

The impact of higher energy costs during the 1970's led to the develop- 
ment of high-efficiency aircraft propulsion systems, most notably the ad- 
vanced turboprop or propfan engine. This engine concept replaced the 
standard fan of the modern turbofan with a high bypass ratio propeller- 
like fan, hence the name propfan. The NASA Lewis Research Center 
(LeRC) development and testing of propfan concepts during the early 
1980's. An unexpected flutter instability occurred for the propfan denoted 
SR-5 and it is reported in Mehmed, Kaza, Lubomski and Kielb [6], This 
propfan was made up of highly swept and twisted flexible blades which 
experienced classical subsonic coupled bending-torsion flutter. Different 
numbers of blades were tested on the rotor and the effect of the increased 
aerodynamic coupling due to higher blade numbers led to a stronger de- 
stabilization of the rotor. 

The occurrence of classical flutter in the SR-5 propfan led to an effort to 
develop more accurate aeroelastic analysis for propfan design. A propfan 
test model was designed by NASA LeRC intentionally to encounter flutter 
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in the wind tunnel. This propfan designated SR3C-X2 (Mehmed and Kaza 
[7]) was a composite blade designed using a specific ply layup to result 
in natural frequencies and mode shapes which would lead to flutter of the 
blade under test conditions. Their work verified the accuracy of the 
aeroelastic analysis methods by demonstrating that a blade could be de- 
signed which would flutter at a specific operating condition. 

1.22 Unsteady Aerodynamic Models 

The underlying physics of an aeroelastic problem requires that the fluid 
pressure response due to the blade vibration be determined: The devel- 
opment of unsteady aerodynamic models for the motion-dependent be- 
havior of cascades can be classified in a variety of ways. The most 
common distinction is made for blades where the influence of airfoil 
shapes can or cannot be neglected. A further complication to the analysis 
is that the unsteady aerodynamic behavior is governed by separate 
equations depending upon if the flow is subsonic, mixed 
subsonic/supersonic (transonic), or fully supersonic. The following dis- 
cussion of unsteady aerodynamic models attempts to cover some of the 
efforts within each of the flow regimes which modern turbomachinery op- 


erates in. 
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One of the first investigations of the inviscid, incompressible unsteady 
aerodynamics of cascades of blades was reported by Whitehead [8,9]. 
His technique was to solve for the unsteady vorticity distribution along flat 
plates which were undergoing rigid translational and rotational harmonic 
oscillations. The airfoil vibration was applied as a boundary condition by 
modifying the velocity upwash along the airfoil surface. The first paper [8] 
was limited to flows over flat plates which caused no steady flow de- 
flection, a truly unloaded flat plate representation of the blades. His later 
work [9] extended the analysis to account for the steady aerodynamic 
loading in the boundary condition gradient terms. The agreement of the 
unsteady aerodynamic results was much better with experimental data 
when the steady loading was included in the unsteady analysis. 

This finding implied that the effect of steady aerodynamic loading was 
important for the unsteady aerodynamics of airfoils, even when modeled 
using simplified, flat plates in incompressible flows. Atassi and Akai [10] 
developed an approach for incompressible flows over airfoils having 
thickness and camber operating within a nonuniform steady flowfield. 
They discovered that airfoil shape and steady aerodynamic loading had 
a very strong influence over the unsteady aerodynamic response of os- 
cillating cascades. 
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Requirements to reduce the noise problem associated with aircraft engine 
fans and compressors presented the need to address compressible flows, 
and the acoustic properties of turbomachinery operating within such 
flows. A development for the unsteady subsonic flow within cascades of 
oscillating flat plates has been reported by Smith [11]. This work was 
primarily concerned with quantifying the acoustical properties of the cas- 
cade by placing special emphasis on the identification of "acoustic reso- 
nance" phenomena for the cascade. The occurrence of acoustic waves 
greatly complicated the compressible unsteady flow problem when com- 
pared to the incompressible flows. 

Namba [12] advanced a similar technique for application to subsonic cas- 
cades of flat plates operating under steady aerodynamic loading due to 
flow incidence. This work concluded that as for the incompressible case, 
the effect of steady aerodynamic loading had a considerable influence 
over the unsteady blade forces. Ni [13] presented a complete unsteady 
aerodynamic analysis for unloaded flat plates for application to both sub- 
sonic and supersonic flows. A thorough discussion of the identification 
of the acoustical characteristics of oscillating airfoils within these flow re- 


gimes was presented. 
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The development of higher tip speed fans and compressors during the 
1970's opened up a unique flow regime for unsteady aerodynamic re- 
search. Operation of turbomachinery at supersonic blade tip rotational 
speeds and subsonic flight speeds results in a special flow problem 
termed "supersonic flow with subsonic leading-edge". This flow condition 
develops when leading edge Mach waves form but do not pass into the 
cascade because the axial flow velocity remains subsonic. A trailing edge 
Mach wave also forms and it impinges upon the adjacent blades in the 
cascade. Considerable effort was devoted to this supersonic flow prob- 
lem because supersonic torsional flutter was becoming more frequent in 
the new high-speed designs. 

Verdon [14] developed a velocity potential method for solving the super- 
sonic flow with subsonic leading-edge problem for a cascade consisting 
of a finite number of oscillating airfoils. The method incorporated an an- 
alytical solution for the potential upstream of the cascade leading edge 
line matched to a finite-difference solution within the cascade region and 
the wake region of the cascade. The Laplace transform solution of 
Kurosaka [15] allowed analysis for infinite cascades of blades, which pre- 
sented considerable improvement in computational time over Verdons 
approach. This formulation resulted in closed-form analytical solutions 
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for the blade unsteady forces which could be solved quickly on the com- 
puter, although the technique was limited to low reduced frequencies. 

Verdon and McCune [16] presented an approach based on a Laplace 
transform solution which had application to the full range of reduced fre- 
quencies. A similar solution using a the Wiener-Hopf technique was de- 
scribed in the work of Adamczyk and Goldstein [17]. This method results 
in a full analytical solution for the unsteady flow problem. Adamczyk and 
Goldstein's model is popular and widely used in the turbomachinery field. 

Miles [18] initiated the unsteady aerodynamic work for flows which had 
supersonic axial flow over the leading edge of the blade in 1956. Lane 
[19] presented a similar analysis which evaluated the integral equations 
and accounted for the Mach wave reflections in the blade-row. Interest in 
supersonic transport in the late 1980's led to renewed interest in the 
analysis of supersonic axial flow. Ramsey and Kielb [20] developed an 
algorithm for solving the formulations of Miles and Lane efficiently on a 
computer. 

Most of the methods discussed to this point have considered the cascade 
to consist of flat plates, or cascades where the steady flow was fully un- 
coupled from the unsteady flow. Ni and Sisto [21] presented an approach 
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which solved the unsteady Euler equations for a cascade of blades using 
a time-marching integration method. The work was limited to flat plate 
airfoils undergoing harmonic motions, but the method was one of the first 
to solve the unsteady Euler equations using a finite-difference solution. 
A key feature of the method was that the mean steady-state flow was fully 
coupled to the unsteady flow in the cascade. A disadvantage of the ap- 
proach was that the time integration resulted in long computer running 
times. 

Another time integration solution for the unsteady Euler equations was 
presented by Huff and Reddy [22], They employed a deforming grid 
method for solving the unsteady Euler equations for cascades of super- 
sonic fan blades modeled using a finite-difference strategy. They pre- 
sented results for a flat plate cascade and compared the solution using 
the unsteady Euler equations and a small-disturbance supersonic theory 
which showed excellent agreement. Results were also included for a 
highly cambered supersonic fan blade undergoing torsional oscillation. 
Disadvantages of this method were excessively long computer running 
times and a limitation which required multiple computational grids for 
small interblade phase angles. 
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Verdon and Caspar [23] developed a solution method for the unsteady 
potential equation by linearizing the unsteady flow about a nonuniform 
steady flowfield. The thickness and camber of the airfoils in the cascade 
was accomodated by solving the unsteady potential as being a small dis- 
turbance about the nonuniform steady potential field. The nonlinear 
steady full potential flow was coupled to the harmonic unsteady potential 
caused by blade vibration. Results showed the effect of airfoil shape and 
flow incidence with comparison to flat plate analysis. An advantage of this 
method was that a direct matrix solution for the unsteady potential was 
obtained compared to iterative solutions commonly used in CFD applica- 
tions. The direct solution scheme results in reduced computational time. 

The transonic unsteady flow problem was investigated by Verdon and 
Caspar [24] in 1984. The linearized unsteady potential flow theory was 
used to model flows containing weak cascade shocks. The shocks were 
assumed to undergo small-amplitude harmonic motion during airfoil os- 
cillation. The unsteady flow was assumed to be a first-order harmonic 
perturbation of the nonlinear steady flow field. Shock fitting was used to 
model the unsteady shock motion. Whitehead [25] presented a solution 
for the linearized unsteady potential equation which was similar to the 
method of Verdon although his approach utilized a finite element solution 
scheme. This technique also used a shock-capturing method for the 
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steady flow, although the results compare very well with those of Refer- 
ence [24], 

1.23 Aeroelastic Analysis Models 

System energy models were commonly used to predict flutter inception in 
turbomachinery prior to 1970. The unsteady aerodynamic loads were 
modeled using isolated airfoil thin-wing theories which were developed 
for aircraft aeroelastic analysis. These isolated airfoil theories were often 
applied at the 3/4 span location of the blade, and an assessment of the 
work that the aerodynamic loads performed could be determined by con- 
sidering the unsteady aerodynamic forces and the blade mode shape. 
Stability was determined based on whether the work of the fluid on the 
blade was either positive (unstable) or negative (stable). An example of 
this type of model is described by Carta [26], His work was principally 
concerned with studying the effect of blade-disk-shroud structural dy- 
namic coupling, although the isolated airfoil theory was incorporated to 
simulate blade-to-blade aerodynamic coupling. 

A need for higher pressure ratios and increased thrust in the 1970's re- 
sulted in turbomachinery having higher tip speeds. Aeroelastic modeling 
utilizing the new unsteady aerodynamic models for turbomachinery more 
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accurately estimated the onset of flutter. Snyder and Commerford [27] 
used an energy method which applied the supersonic/subsonic leading 
edge model of Verdon [14] to predict flutter of high-speed fan flutter. In 
addition, they also performed an experiment in a linear cascade of seven 
oscillating blades to measure when flutter would occur for ranges in Mach 
number and frequency. The agreement between the supersonic unsteady 
aerodynamic model and the measured flutter conditions was very good. 
It was noted that the cascade effects of the adjacent blades caused the 
system to be less stable than the isolated airfoil theory predicted, and 
they concluded that cascade unsteady aerodynamic theory must be used 
for the unique unsteady flows of turbomachinery. 

A description of the flutter prediction model used by Pratt & Whitney Air- 
craft in 1975 is outlined in the work of Mikolajczak, et. al. [28], This work 
discovered that the empirical correlations and design rules formerly used 
were inadequate for flutter analysis of high-speed flows. The aeroelastic 
model was an energy model where the sum of the calculated aerodynamic 
and structural damping was used as a measure of stability. The model 
incorporated unsteady aerodynamic theories of Whitehead [8], Smith [11], 
and Verdon [14] for estimation of unsteady aerodynamic forces in the 
incompressible, subsonic, and supersonic flow regimes. This work em- 
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phasized the importance of accurate blade vibrational analysis, blade-disk 
structural coupling, and unsteady aerodynamic forces in flutter analysis. 

The aeroelastic models up to 1980 were primarily based on energy meth- 
ods, where the blades were assumed to vibrate in a single degree-of- 
freedom or in a single vibrational mode. A new trend in flutter analysis 
was to include multiple degrees-of-freedom for the blade and to investi- 
gate the effect of coupling on blade flutter. Bendiksen and Friedmann [29] 
presented an aeroelastic model which incorporated two degrees-of- 
freedom (DOF) for the blade. This approach utilized the incompressible 
cascade theory of Whitehead [8], A lumped parameter model of a "typical 
section" along the span of the blade was used which had one bending and 
one torsional degree-of-freedom. The results imply that bending-torsion 
coupling can have a pronounced effect on cascade flutter, and that accu- 
rate aeroelastic analysis must account for multiple degrees-of-freedom. 

The lumped parameter two DOF model became the most commonly used 
approach for aeroelastic analyses during the 1980's, and serves as the 
workhorse of modern turbomachinery flutter analysis. The combination 
of the unsteady aerodynamic theories of Smith [11], Adamcyzk and 
Goldstein [17], and Ramsey and Kielb [20] combined with the coupled 
bending-torsion two DOF blade description has resulted in flutter predic- 
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tion models which are relevant for all of the flow regimes of modern 
turbomachinery. 

All of the aeroelastic models discussed assume that each blade within a 
rotor is identical, having exactly the same natural frequencies and mode 
shapes. This system is termed a "tuned" system. Research to determine 
the aeroelastic behavior of systems which have small statistical variations 
in dynamic properties of the adjacent blades on the rotor has begun. 
Such a rotor is referred to as "mistuned" because the dynamic behavior 
of the rotor is no longer uniform from blade-to-blade. A study of mistuned 
systems ([30] to [33]) has concluded that the tuned system represents the 
most unstable configuration for a rotor. Mistuning has a stabilizing effect 
on the flutter stability of rotors, primarily due to a break-up of the cyclicity 
of the rotor. This was an important conclusion because the analysis of 
mistuned rotors was computationally expensive, and knowledge that the 
tuned analysis is conservative simplifies the turbomachinery designer s 
task. 

The occurrence of flutter in the SR-5 advanced propfan indicated that a 
flutter problem existed which could not be predicted using the available 
aeroelastic analysis methods. The propfan was a thin, swept, flexible 
propeller blade which experienced large-amplitude self-excited vibration 
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during performance testing at NASA LeRC [6], The complicated 
vibrational behavior of the blade could not be properly predicted using a 
simple typical section two DOF aeroelastic model. 

Nonlinear large-displacement beam structural dynamic models were in- 
troduced in order to incorporate the three-dimensional properties of the 
blade. Kaza and Kielb [34], [35] developed an analysis which modeled 
the blades as cantilevered beams having varying properties along the 
span. The model also accounted for the effects of centrifugal stiffening of 
the beam. The unsteady aerodynamic forces were evaluated at several 
radial positions along the span of the beam and the forces were numer- 
ically integrated along the span to arrive at an unsteady aerodynamic 
generalized force. The unsteady aerodynamic models used were those 
of Smith [11] for subsonic flows and Adamczyk and Goldstein [17] for 
supersonic flows. 

The beam models of the propfans resulted in better analytical prediction 
of the experimental flutter conditions. But the wind tunnel propfan models 
underwent complicated coupled vibrational modes, which could not be 
fully simulated using the beam representation. The finite element method 
was chosen to model the propfan blades so that the simulation of the 
flexible blade dynamics could be accomplished more accurately. 
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Turnberg [36] first applied the finite element method to the analysis of 
propfan flutter in 1983. He used modal analysis to base the aeroelastic 
response on the in-vacuum free vibration modes of the blade. The sub- 
sonic theory of Smith [11] was used to determine the unsteady aerodyna- 
mic forces acting at several "strips" along the span of the blade. The 
comparison with experimental results was very good. A disadvantage of 
the method was that the aerodynamic model assumed that the flow on 
each strip was two-dimensional, even though experimental flow visual- 
ization found that large spanwise flows occurred for the propfan blades. 

An extension of the finite element based modal aeroelastic analysis was 
developed by Kaza, et. al. [37] by incorporating a three-dimensional un- 
steady aerodynamic theory. The aerodynamic model used a lifting sur- 
face theory to model the blade using oscillating doublets which were 
placed on panels along the three-dimensional surface of the blade. This 
approach showed excellent agreement with the wind tunnel experiments 
for flutter prediction over a wide range of speeds. The primary limitation 
of the lifting surface aerodynamic model was that the airfoil was modeled 


as a zero-thickness camber line with zero flow incidence. 


- 22 - 


1.3 Objective of Current Work 

The purpose of this work is to present the development and application 
of an aeroelastic stability analysis which is suitable for all types of axial- 
flow turbomachinery. A review of the related previous work has discov- 
ered most of the aeroelastic analyses developed to date have been 
restricted to two degree-of-freedom lumped parameter models or flat 
plate small-disturbance unsteady aerodynamic models which disregard 
the effects of steady aerodynamic loading. The unsteady aerodynamic 
models which account for steady aerodynamic loading often require con- 
siderable computer running times which prevent their use within pro- 
duction aeroelastic codes. Based on the results of the literature review, 
the present work was begun with the objectives to: 

• Develop an aeroelastic stability analysis suitable for general axial-flow 
turbomachinery blade designs having varying degrees of airfoil 
camber and thickness and operating in compressible flows. 

• Utilize two-dimensional steady full potential and linearized unsteady 
potential flow theory to calculate the motion-dependent aerodynamic 
loads acting on the blade during vibration. 
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• Account for three-dimensional aeroelastic effects by calculating the 
unsteady aerodynamic loads on two-dimensional strips which are 
stacked along the span of the blade from airfoil hub to tip. 

• Use modal analysis to couple the in-vacuum natural vibrational modes 
of the blade to the motion-dependent unsteady aerodynamic loads 
within the frequency domain. 

• Demonstrate capability by applying method to a high-energy subsonic 
turbine rotor blade which is suspected to experience aeroelastic 
problems. 

The implementation of this analytical formulation was performed at NASA 
LeRC in support of their turbomachinery aeroelasticity research. NASA 
LeRC is developing an aeroelastic analysis for turbomachinery called the 
Forced REsponse Prediction System (FREPS) which is applicable to flutter 
and forced response prediction for fans, compressors, and turbine blades. 



CHAPTER 2 - ANALYTICAL DEVELOPMENT 


2.1 Aerodynamic Analysis 

The fluid flow within the turbomachinery blading is modeled along stacked 
two-dimensional streamsurface "strips" over the blade from the airfoil hub 
to tip. This model is a simplification of the three-dimensional flows which 
occur in turbomachinery blades. However, it is fair to assume that the 
streamsurfaces for a high hub/tip ratio turbomachine are almost two- 
dimensional, except close to the hub and the tip of the blade where sec- 
ondary flows dominate. 

The linearized potential theory of Verdon and Caspar [24] is used to 
model the unsteady flow within the blades. This linearized potential ap- 
proach assumes that the unsteady flow resulting from airfoil oscillation is 
a small perturbation of the nonlinear steady-state (or "mean") flow. It is 
assumed that the primary influence of blade thickness and camber on the 
unsteady flow is due to the motion of the airfoil through the nonuniform 
steady flowfield. The steady full potential flowfield is evaluated prior to 
solution of the unsteady potential flow in order to account for the coupling 
of the steady and unsteady flows. A brief description of the governing 
equations and solution for both the steady full potential flow and the 
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linearized unsteady potential flow follows. A detailed review of the steady 
and unsteady aerodynamic formulation used may be found in Verdon [38]. 

2.11 Steady Full Potential Flow 


The steady compressible flow along a two-dimensional surface is calcu- 
lated using the full potential method of Caspar, et. al. [39]'. A finite-area 
technique is utilized to solve the continuity equation for a single blade 
passage. For compressible, inviscid, isentropic, and irrotational steady 
flow of a perfect gas, the continuity equation can be solved uniquely for 
the steady flow. 


V-(pl/) = 0 (1) 

In terms of the velocity potential, equation 1 yields the steady full potential 
equation (equation 2). 


V-(pV<I>) = 0 


( 2 ) 


The flow variables p and V are nondimensionalized by their values at the 
inlet of the passage denoted as state 1. This nondimensionalized 
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equation, along with the steady Bernoulli equation and the ideal gas 
equation of state results in an expression for the local density. 


(y-i) 


/If, (VO VO 


-] 


_J 

y- 1 


( 3 ) 


The inlet Mach number at state 1 is M, and the specific heat ratio of the 
gas is y. The steady potential flow within the cascade may be solved us- 
ing the continuity equation of eq (2) and the local density description of 
eq (3). The solution is also solved for weak shock waves for which the 
flow is assumed to remain isentropic and irrotational. The steady flow 
analysis incorporates an artificial viscosity scheme in supersonic regions 
to stabilize the calculations (by Caspar [40]). 


Figure 1 includes a figure showing a two-dimensional representation of 
the cascade with the inlet and exit flow velocity vectors indicated. A dia- 
gram of the cascade geometry, including the stagger angle, airfoil chord, 
cascade gap and the leading and trailing edges is included on figure 2. 
The regions where the boundary conditions for the full potential equation 
are applied are shown on the blade-to-blade cascade of figure 3. 

The boundary conditions for the mean potential flow require the airfoil 
surface flow tangency condition (equation 4). 
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VO n = 0 along the airfoil surface (4) 

An additional boundary condition is applied at the periodic boundaries 
as 


lower periodic upper periodic . 

which states that the flow is continuous across the periodic boundaries. 
The boundary condition at the upstream boundary 1 requires specification 
of the inlet Mach number /Vf, and the inlet flow angle /?,. 

The airfoil is assumed to have a sharp trailing edge, such that the Kutta 
condition may be satisfied. The addition of a wedged trailing edge also 
simulates the viscous interaction of the upper and lower surface flows 
which meet at a blunt trailing edge. This Kutta condition allows the 
downstream flow angle /? 2 to be prescribed as the bisector of the blade 
trailing edge metal angle. The overall mass conservation can then be 
solved for the downstream exit Mach number M 2 . 

The finite-area solution of equation (2) requires the use of a spatial mesh 
to discretize the flow passage. The approach of Reference [40] uses two 
spatial meshes for the mean flow solution. An H-type of computational 
mesh is used to capture the overall flow details of the passage. A C-type 
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of computational mesh is also used to resolve the finer flow details 
around the leading edge (LE) stagnation point and around shocks. Ex- 
amples of these types of meshes are shown on figure 4. The mass con- 
servation equation is applied to the cells of the finite-area mesh for the 
solution of equation (2). A successive line over-relaxation scheme is used 
to iteratively solve the nonlinear steady potential equation. 

2.12 Linearized Unsteady Potential Flow 

The unsteady flow, i.e. that flow caused by airfoil oscillation, is calculated 
using the linearized potential approach of Verdon and Caspar [24]. The 
unsteady flow is modeled as a small-amplitude potential perturbation 
from the steady potential flow. The unsteady potential is a scalar repre- 
sentation of the velocity which requires that the unsteady flow remain 
irrotational and isentropic. The airfoils in the cascade are assumed to 
vibrate with small-amplitude harmonic motion and a constant phase angle 
between adjacent blades. 

These assumptions permit the unsteady equations of motion for the flow 
to be reduced to a single partial differential equation which is solved for 
the harmonic unsteady potential. The unsteady potential equation is 
shown in nonconservative form in equation 6. 
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In this equation <t> is the steady velocity potential and <f> is the unsteady 

velocity potential. The unsteady substantial derivative operator is 

= ico() + VO f). Derivation of equation (6) from the unsteady Euler 
Dt v/ 

equations is shown in Appendix A. 

The boundary conditions for the unsteady potential problem require 
specification of conditions along the upstream and downstream, vibrating 
airfoil surfaces, and the cascade periodic boundaries (figure 5). The sur- 
face flow tangency condition requires that the unsteady flow remain 
bounded by the moving airfoil surfaces. The airfoil surface may undergo 
any general deformation mode, as indicated in figure 6. The surface 
boundary condition along the airfoil is 

V0n= -|j- + V T -^--(fV)Voj-n 


( 7 ) 
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where n is a unit vector normal to the airfoil surface, r is the airfoil dis- 
placement vector, V, is the velocity tangent to the airfoil surface, and r is 
the unit vector tangent to the airfoil surface. 

The first term of equation (7) defines the velocity of the airfoil normal to 
its surface during oscillation. The second term represents the rotation of 
the airfoil through the steady potential field. The last term extrapolates 
the surface boundary conditions along the surface of the oscillating airfoil 
to the mean location of the airfoil in the steady reference frame. 

Expressing the surface boundary conditions of the oscillating airfoil at the 
steady airfoil location permits solution on a stationary, nondeforming 
computational grid. A Taylor series expansion of the moving airfoil lo- 
cation relative to the steady airfoil location is applied to extrapolate the 
boundary conditions to the steady airfoil surfaces. First-order (linear) 
terms are retained in this series expansion (Reference [38]) which results 
in the last term of equation (7). 


Substitution of the harmonic time dependence into equation (7) results in 
the surface boundary condition of equation (8). 


V<f)-n = [/'wr 4 - (VO-T)(r-V)r — (r V)V<t>J n 


( 8 ) 
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The periodic boundary condition for the unsteady potential is prescribed 
as a constant phase angle difference in time for the adjacent blade in the 
cascade. Since the airfoil is assumed to undergo harmonic motion, the 
solution for flow passages adjacent to the computational passage simply 
vary by a constant harmonic "interblade phase angle a (see figure 5). 
The adjacent blade's harmonic motion in terms of the interblade phase 
angle is given by equation (9). 

7 = (9) 

The j represents the blade number on the rotor. All of the flow variables 
are assumed to vary between adjacent flow passages by this interblade 
phase angle relationship. The upstream and downstream boundary con- 
ditions are determined from analytical expressions for the far-field poten- 
tial fluctuations due to blade oscillation. 

The equation for the unsteady potential is solved using a weighted least- 
squares finite difference discretization. Differencing weights are deter- 
mined from the coefficient terms of equation (6). The blade passage is 
discretized using the same type of global and local computational meshes 
as described for the steady potential computation (figure 4). Solution for 
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the unsteady potential is obtained by a direct matrix solution of equation 

( 6 ). 

The unsteady pressure P is found from the unsteady expansion 


P = p + (r V)P (10) 

where p represents the harmonic pressure and the second term repres- 
ents the first-order Taylor series expansion for the pressure variation due 
to the airfoil motion through the steady potential and pressure fields. 

The harmonic pressure is calculated from the unsteady potential as 

_ D4> 


p = — p(/CO0 - |- (11) 

The unsteady forces and moments acting on the airfoil are determined by 
integrating the harmonic pressure of equation (10) over the airfoil surface. 

The unsteady flow problem described in this section is determined for a 
prescribed set of steady flow and unsteady flow parameters. The steady 
flow parameters for a particular airfoil geometry and working fluid are (a) 
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the inlet Mach number /W, and (b) the inlet flow angle /?,, where the exit 
Mach number and angle are determined from the trailing edge Kutta 
condition. The unsteady flow parameters are (a) the blade displacement 
vector r, (b) the interblade phase angle a, and (c) the vibrational reduced 
frequency k = — — . 

The nondimensional reduced frequency k is a measure of the unstead- 
iness of the flow. In general, the reduced frequency can be viewed as the 
ratio of the time a fluid particle takes to pass over the surface of the blade, 
i.e. proportional to — , compared to the time it takes to complete one cycle 
of blade vibration Thus, when the reduced frequency is low, i.e. less 
than or equal to 0.5, the flow is varying in approximately the same fashion 
as the excitation. This situation can lead to large unsteady aerodynamic 
forces. Likewise when the reduced frequency becomes high, i.e. larger 
than 0.5, the excitation is much higher than the fluid particle velocity and 
the flow has little time to react to the excitation. Flows at higher reduced 
frequencies usually result in relatively low unsteady aerodynamic forces. 

2.2 Structural Dynamic Analysis 

Dynamic analysis of blades of general shape is most efficiently accom- 
plished using the finite element method (FEM). The FEM permits a 
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straightforward method for estimating the dynamics of blades having 
complex geometries and being made of advanced anisotropic materials. 

The dynamic equations of motion for the finite element representation of 
the blade are shown below 

+ [C]{u} + [/OW = {/>«} + {/=„(!)} (12) 

where the n global degrees-of-freedom (DOF) displacements are {u}. The 
structure mass, damping, and stiffness are represented by the matrices 
C/VfJ,[C], and f/Cj, respectively. The forcing terms on the right-hand 
side of equation (12) represent the external forces due to aerodynamic 
{F A (t)} and mechanical {f M (f)} sources. 

The undamped homogeneous form of equation (12), i.e. neglecting forcing 
terms and damping, represents the free vibration problem of the blade. 
Solution of the free vibration problem requires the solution of the real 
eigenvalue problem 


a) 2 lM]{u} = [K]{u} (13) 

which results in the set of n eigenvalues co n and eigenvectors {</>„}• These 
eigenvalues represent the undamped natural frequencies for the blade 
and the undamped natural modes, or mode shapes. 
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The dynamic motion of the blade, in the presence of forcing functions and 
damping, is assumed to be represented as a linear combination of con- 
tributions from each of a reduced set of the free-vibration natural modes. 
This modal expansion using m modes for the n DOF system assumes that 
the dynamic displacements can be written as 

{u(t)} = W{q(t)} (14) 

where [O] is the n x m modal matrix having the undamped eigenvectors 
{</>„} arranged columnwise. The vector (q(f)} is the time-dependent modal 
coordinate vector. The modal coordinates represent a measure of each 
of the natural vibrational mode's contribution to the dynamic displace- 
ments. The choice of which m modes to retain for the dynamic analysis 
depends upon the nature of the expected forcing function. The contrib- 
ution of modes having high natural frequencies relative to the forcing fre- 
quency will be low. In general, only those modes having frequencies 
close to the expected forcing frequency or vibrating frequency need to be 
retained. 

Substitution of equation (14) into the original dynamic equations of motion 


equation (12), and premultiplying by the transpose of the modal matrix 


- 36 - 


transforms the dynamic equation from physical space to modal 
space. The modal equation of motion becomes 

F M G ]{q} + [C G ]{q} + £K G 3{q} = [4>]V/i(')} + {F M (t)}) (15) 

where this equation represents m equations, so that the dynamic problem 
has been reduced from a system of n equations. 

The following definitions are used for the transformed matrices 

[C G ] = [<D] r [C][®] (16) 

These matrices are known as the generalized mass, generalized damping, 
and generalized stiffness matrices and they are of size m x m. The gen- 
eralized mass [M a ] and generalized stiffness (Xd matrices are diagonal 
due to the orthogonality of the eigenvectors to the physical mass [Mj and 
stiffness [X] matrices. 

It is difficult to accurately determine the generalized damping matrix 
[C G ] because of the complicated nature of damping in turbomachines. 
Damping for a blade is strongly dependent upon material damping, 
damping in connections, such as contact stresses in turbine blade roots, 
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and structural interfaces such as shroud gaps. These types of con- 
nections are often a nonlinear function of displacements, and a consider- 
able amount of research has been devoted to damping models with little 
definite success. 

A common assumption for modeling generalized damping is to assume 
equivalent modal damping, which is similar to single DOF damping, ap- 
plied to each individual vibrational mode of the structure. The concept 
of modal damping results in the following diagonal generalized damping 
matrix 


[C G ] = (17) 

where the modal damping ratios ( corresponding to each individual mode 
are along the diagonal. 

The force expression on the right-hand side of equation (15) is referred to 
as the generalized force vector {0(f)} and it is defined in equation (15). 

{0(0} = [<p] «f 4 (0! + {wm (is) 


Using the definitions 
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{<?„(')} = m W)> 

my = w Ww 

simplifies the total generalized force as the sum of the aerodynamic gen- 
eralized forces (Q*(f)} and the mechanical generalized forces (Q M (f)} 

The final form of the modal equations of motion, using the above defi- 
nitions, becomes 


+ CC G D{qr} + £K G 3{q} — {Q(f)} (20) 

which represents the m modal equations for the system. Note that this 
equation, with the assumption of modal damping, represents a system of 
m uncoupled, second-order ordinary differential equations. This is a sig- 
nificant simplification from the system of n fully coupled equations as 
given by equation (12). The loss of accuracy due to the modal truncation 
can be reduced by including more modal coordinates, and basis mode 
shapes, within the analysis. 

2.3 Aeroelastic Model 


The rotor for the present study is assumed to vibrate as a "tuned" rotor, 
where every blade on the disk has the same natural frequencies and 
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mode shapes. This assumption is obviously inaccurate, because every 
blade will have small differences due to the manufacturing process, in- 
stallation effects, etc. But a tuned rotor always results in the most con- 
servative estimate of flutter, as was discussed in section 1.23, and 
additionally significantly simplifies the computational effort for aeroelastic 
analysis. 

The tuned rotor implies that the rotor blades have equal amplitudes of 
motions while having a constant phase lag between adjacent blades. This 
assumption results in a travelling wave which passes along the periphery 
of the rotor with a certain frequency and wavelength which corresponds 
to the constant interblade phase angle a. A tuned rotor may vibrate with 
any of N possible interblade phase angles where N is the number of 
blades on the rotor. 


The interblade phase angles are limited to the discrete values 


2rc(/ - 1) 

N 


for j = 1 to N. The vibration of the rotor in the <j‘ h interblade phase angle 
refers to the j ,h travelling wave of the rotor. 
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Incipient flutter is a linear problem since the vibration up until it becomes 
unstable is small-amplitude and occurs at a definite individual frequency. 
Flutter prediction is concerned with identifying the conditions where the 
vibrations just begin to go unstable. 

Since the vibrations occur at a single frequency, the modal coordinate 
vector can be assumed to vary harmonically in time. This is represented 
as 


M0} = {9 (22) 

where the vibrational frequency, i.e. the flutter frequency to,, is assumed 
to be close to a natural frequency of the blade. Such an assumption im- 
plies that the dynamic displacements of the aeroelastic system {u(f)} are 
harmonic, and consist of a linear combination of harmonic vibrations of 
each undamped mode which corresponds to the modal coordinates {q}. 

The unsteady aerodynamic forces on the blade are applied to the blade 
as the external aerodynamic forces (F^(f)}. When the mechanical forces 
are neglected, the generalized forces reduce to only the aerodynamic 
generalized forces (equation 23). 

{ 0 ( 0 } = m T {F A (t)} 


(23) 
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The external aerodynamic forces may be classified as either (a) motion- 
dependent aerodynamic forces or (b) motion-independent aerodynamic 
forces. The motion-dependent forces are caused by the response of the 
fluid to the underlying blade motion. The pressure response of the fluid 
lags the motion of the blade, so that the motion-dependent forces are re- 
presented as complex having in-phase (real) and out-of-phase (imaginary) 
parts. The motion-independent forces result from aerodynamic 
excitations such as incoming pressure waves or velocity wakes, and these 
forces are primarily of importance in determining the forced response of 
aeroelastic systems. For linear flutter analysis, they play no role in the 
prediction of stability and will be neglected. 

The motion-dependent unsteady aerodynamic forces are modeled as 
complex nodal loads which act at the finite element physical DOF. Since 
linear unsteady aerodynamic theory was assumed the full unsteady forces 
can be represented as a summation of the forces caused by motion at 
each of the natural modes of vibration of the blade. This assumption is 
similar to the expansion of the physical DOF using modal analysis. 

Following this assumption, the motion-dependent unsteady aerodynamic 
forces are modeled as 


(MOl = «, 2 [R]{q(0! 


( 24 ) 
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with [ft] being the n x m unsteady aerodynamic modal forces having 
each individual force vector {ft} arranged columnwise, as in the modal 
vector [d>]. The force vector {ft}, is the unsteady aerodynamic forces 
acting at each blade DOF due to vibration of the blade in the j ,h mode. The 
coefficient co, 1 2 3 4 results because of the scaling of the unsteady pressures by 
the inlet steady dynamic pressure term pV, 2 . 

The modal unsteady forces are strongly dependent upon the assumed 
flutter frequency oj t , interblade phase angle a, inlet Mach number /Vf, and 
the rotor speed. The calculation of the modal aerodynamic forces {ft}, is 
accomplished in the following procedure: 

1. Interpolate the blade mode shape {<£}, from the finite element model 
to the unsteady computational mesh. 

2. Solve the linearized unsteady potential equation using the mode 
shape {<f>}, as the airfoil displacement vector r. 

3. Calculate the resulting unsteady pressure p along the airfoil surface. 

4. Integrate the surface unsteady pressure p along the airfoil surface to 
result in unsteady concentrated forces {ft}, which act at the finite ele- 


ment nodal DOF. 
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The resulting force vector due to the j th airfoil mode shape, is the modal 


unsteady force vector {R}, 


Substitution of the motion-dependent force expression of equation (24) 
into the generalized force of equation (23) results in the following 

(Q(f)} =a), 2 [(D] r [«]{q(f)} 

or 

{O(0> = CO, 2 [/I]{<»(0> (25) 

and the modal aerodynamic matrix is 

[4] = Oft/?] (26) 


which is of size m x m. 

A physical interpretation of [A] is that the A u element represents the 
aerodynamic generalized force in the i lh mode caused by motion of the 
blade in the j' h vibrational mode. This matrix is in general nonsymmetric, 
and has complex elements which include the magnitude and phase of the 
generalized unsteady aerodynamic forces. The matrix also requires har- 
monic time variation because the unsteady aerodynamic model used is 
only applicable for harmonic blade motions. 
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Definition of the modal aerodynamic matrix permits substitution of eqs 
(25) and (26) into the equation of motion equation (20) to result in 

CM g J{c/} + rC G J{q} 4- = <*>f 2 LA]{q}. (27) 

The harmonic motion assumption could now be applied to reduce this 
system to an eigenvalue problem. Unfortunately, the presence of damp- 
ing in this equation prevents a straightforward eigensolution of the modal 
equations because they are in quadratic form. Eigenvalues of a damped 
system occur in complex conjugate pairs, resulting in 2m eigenvalues, 
whereas the modal equation (27) contains only m equations. 

A state vector transformation is used to express equation (27) in state 
space, which results in a system of 2m equations for the 2m eigenvalues. 
This transformation (Reference [41]) involves defining the auxiliary vari- 
ables 


(0 = 



(28) 


which results in the state space representation of the system of equations 


= [«']«) 


( 29 ) 
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with the modified partitioned matrices defined as 


[M‘] = 


[o] (£ m g 3 + C/tK)]) 

([TM G J + C^Oty)]) 


[*’] = 


(OU + OH)]) Co] 
[o] - Dk'gJ 


(30) 


The state-space variables {£} are assumed to have harmonic variation in 
time, as 

{{}-{«}«" ( 31 > 


so that the transformed modal equation of equation (29) becomes the fol- 
lowing aeroelastic eigenvalue problem 


= [*]{{} 


(32) 


which is a 2m x 2m complex eigenvalue problem. 


The set of complex eigenvalues from the solution of equation (32) are ex- 


pressed as 
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). = n + v; (33) 

where the real part (n) of the eigenvalues represents a measure of the 
system damping, both mechanical and aerodynamic, and the imaginary 
part (v) of the eigenvalues represents the damped natural frequency of the 
system. 

The '2m' eigenvalues are used to assess the stability of the system for the 
vibration at the frequency c a,. The stability of the system is governed by 
the real part of the eigenvalues. When the real part n is greater than or 
equal to zero, the system will be unstable, with the vibration amplitude 
growing exponentially in time. Flutter is said to occur when n = 0, a 
neutrally stable condition, or when n > 0, an unconditionally unstable 
condition. 

The eigenvalue problem represented by equation (32) requires an itera- 
tive solution because of the dependence of the modal aerodynamic matrix 
upon the assumed frequency c o f . Solution of this problem, for a specific 
rotor speed and Mach number, requires calculation of the modal aero- 
dynamic matrix and damped eigenvalues for an assumed frequency co f . 
The calculated eigenvalue of interest from equation (32) must have the 
frequency portion v equal to the assumed flutter frequency co, because of 



- 47 - 

the restrictions of the unsteady aerodynamic theory. If the calculated 
eigenvalues are not equal to the assumed flutter frequency, the new as- 
sumed flutter frequency is set to the calculated damped natural frequency 
co,= v. The modal aerodynamic matrix is recalculated and the 
eigensolution continued in this fashion until convergence is achieved. 



CHAPTER 3 - APPLICATION OF METHOD 


3.1 Computer Program 

An outline of the computer implementation of the aeroelastic model from 
Chapter 2 is diagrammed in figure 7. The procedure involves reading in 
the blade geometric information, i.e. finite element model, and also read- 
ing the free-vibration eigenvalues and eigenvectors. The free-vibration 
analysis is performed using either the MSC/NASTRAN or the MARC 
general-purpose finite element packages. The aeroelastic strip definitions 
are input after the blade geometry and modal information. Each strip is 
defined by specifying a locus of finite element nodal points along the 
airfoil surface. For each strip the inlet velocity triangle for the flow is 
specified. Additionally, the aeroelastic parameters (flutter frequency co f , 
interblade phase angles a) are input. 

Along each strip, the computational grids are generated for the calcu- 
lation of the steady and unsteady potential flows. The method currently 
uses two different sets of grids for the steady and the unsteady flow 
problems. The steady flow solution for each strip is calculated and stored 
in database files on the computer. Similar information is initialized for the 
unsteady flow on each strip and stored in the database files. 
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FREPS begins the aeroelastic analysis by initializing the supplied data- 
base files which contain the steady potential solution and the computa- 
tional grid information. The modal aerodynamic matrix is calculated 
using the procedure outlined in figure 8. The airfoil two-dimensional 
mode shapes are determined from the finite element eigenvectors for 
each strip. The mode shapes are interpolated from the finite element ref- 
erence frame onto the unsteady aerodynamic computational grid. The 
unsteady flow due to the enforced airfoil mode shape is solved along each 
strip by the method of section 2.12. 

Numerical integration of the unsteady pressures along the airfoil surface 
is used to calculate the modal unsteady aerodynamic force vector {R}, 
due to the j ,h natural mode. The j' h column of the modal aerodynamic 
matrix [/t] is evaluated in a loop over all the modes as A„ = {4>} T ,{R},. 
This procedure is followed for all of the blade natural modes to completely 
determine the matrix. 

3.2 SSME HPOTP Turbine Description 


The method outlined above was used for the aeroelastic stability analysis 
of the Space Shuttle Main Engine (SSME) turbopump power turbine. The 
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SSME is a hydrogen-fueled liquid rocket engine which generates 512,000 
Ibf of thrust at full power. Each rocket engine is fed by four turbopumps, 
two high pressure and two low pressure booster pumps. The high pres- 
sure turbopumps are driven by axial-flow gas turbines which are powered 
by hydrogen-rich steam generated in individual preburners. 

The High Pressure Fuel TurboPump (HPFTP) supplies the liquid hydrogen 
propellant to the main combustion chamber. The power turbine for the 
HPFTP generates 75,000 horsepower at a rotor speed of 36,000 RPM. 
Each turbine blade of this rotor transmits over 700 horsepower. The High 
Pressure Oxidizer TurboPump (HPOTP) supplies the oxidizer to the com- 
bustion chamber. The HPOTP gas turbine generates 30,000 horsepower 
at 28,000 RPM with each blade transferring approximately 300 horse- 
power. 

Both the HPFTP and the smaller HPOTP operate in a severe environment 
at pressures of up to 5000 psi and temperatures close to 1600 deg F. 
Turbine blade cracking has been a continual problem for both of the high 
pressure turbopumps. A material replacement of single-crystal alloys for 
the HPFTP has effectively reduced the blade cracking problem. 
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A cross-sectional view of the HPOTP turbopump is provided on figure 9 
which shows the pump impellers, shaft, gas turbine, preburner, and the 
turbine blade coolant jet ring. The first stage turbine blade has experi- 
enced frequent cracking in the shank region of the blade. The cracking is 
believed to be due to a resonant excitation of the second vibrational mode 
of the blade with cooling jets which direct a steam mixture to the blade 
shanks. There are 19 equally spaced cooling jets distributed circumfer- 
entially at a radius on the shank just above the top fir-tree attachment 
lobe. The blade dimensions are approximately 1.4 in. (3.56 cm) tall, tip 
chord of 0.66 in. (1.68 cm) and an airfoil span of 0.5 in. (1.27 cm) with a tip 
diameter of 10.8 in. (27.43 cm). 

The HPOTP first stage turbine has blade-to-blade friction dampers in- 
stalled to provide additional mechanical damping to keep vibrational am- 
plitudes low. The original friction dampers were a one-piece design which 
did not provide sufficient damping to prevent fatigue cracking early in the 
development program. A redesign of the dampers in 1985 to a led to a 
reduction in cracking and longer service life, but the desired design life 
of 55 launches has yet to be realized. 

NASA LeRC's involvement in the SSME HPOTP blade cracking problem 
was to assist in identifying methods for determining the level of mechan- 
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ical and aerodynamic damping which exist on the HPOTP turbine blade 
and to estimate the vibratory stress levels caused by the cooling jet and 
the blade-row aerodynamic interaction forcing functions. The aerodyna- 
mic damping for these blades is essentially due to the motion-dependent 
unsteady aerodynamic response of the fluid as the blade undergoes vari- 
ous modes of vibration. This report presents the aeroelastic stability 
analysis for the HPOTP turbine blade, including estimates of the level of 
aerodynamic damping, and a determination of the stability of the motion. 

3.3 Aeroelastic Model 


The aeroelastic model for the HPOTP first stage turbine blade consists of 
(a) the finite element model of the blade and (b) the aerodynamic strip 
definitions along the airfoil of the blade. The finite element model used 
for this work was provided by the the SSME contractor, Rocketdyne Divi- 
sion of Rockwell International Corp. The original finite element model was 
from the ANSYS general-purpose finite element program, and the model 
was converted into MSC/NASTRAN form for analysis at NASA LeRC. 

The finite element model is three-dimensional, consisting of 10,014 nodal 
points and 7758 solid hexahedron elements as shown on figure 10. The 
turbine blade has extensions for the blade-to-blade friction dampers lo- 
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cated on the sides of the blade, between the top firtree lobe and the blade 
platform. This blade has a hollow core, to reduce the blade mass and 
thermal inertia and a tip seal to prevent spanwise flows over the tip of the 
blade. 

Six aeroelastic strips, which represent approximate streamsurfaces, are 
defined along the span of the airfoil, from the blade platform to the tip 
shroud. These strips coincide with sections of the finite element model 
which have constant radii for the nodal points. The strips are indicated 
on figure 11, and a table of the cascade properties for each of the strips 
is included as Table I. 


A plot showing the airfoil cross-section of the HPOTP turbine is included 
on figure 12, along with a table of airfoil coodinates. This turbine blade 
has blunt, rounded leading and trailing edges. The maximum airfoil 
thickness is approximately 30 percent of the airfoil chord length. The 
turbine turns the relative flow into the blade through up to 130 degrees. 
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3.4 Structural Dynamic Results 

The MSC/NASTRAN finite element program was used to calculate the in- 
vacuum blade natural frequencies and modes. Springs were applied at 
the eight load bearing surfaces of the fir-tree lobes to simulate disk flexi- 
bility. The spring constants for these springs were determined by match- 
ing the first three calculated natural frequencies with those measured by 
Rocketdyne during high-speed rotating rig dynamic testing. 

The effect of rotational speed was incorporated within the analysis by 
adding centrifugal loading and a temperature variation along the span of 
the blade. The temperature distribution for several operating speeds was 
estimated based upon flow path gas temperature measurements obtained 
during actual engine operation. The thermal variation of the blade mate- 
rial properties, for the directionally solidified alloy MAR-M-246, was 
modeled using a tabular description of the orthotropic material constants 
versus temperature. 

The natural frequencies were calculated by performing a geometric non- 
linear static analysis (MSC/NASTRAN solution sequence 64) to simulate 
blade rotation and thermal loading. The differential stiffness matrix from 
the static solution was stored in a database file. The normal modes anal- 
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ysis (MSC/NASTRAN solution sequence 63) was performed by including 
the differential stiffness matrix within the linear global stiffness matrix. 
This method of accounting for the effect of centrifugal and other combined 
loadings on the normal modes of rotating structures is described by 
Lawrence, et. al. [42]. 


The first four natural frequencies of the blade were calculated at pump 
rotational speeds of 0, 19,500 and 28,000 RPM, speeds representative of 
typical power levels for an actual shuttle launch. A Campbell diagram is 
included as figure 13 which shows the variation of the natural frequencies 
with increasing pump speed. This blade shows only a very weak de- 
pendence of frequencies on the rotational speed. The effect of rotational 
speed on the frequency is due to centrifugal stiffening of the blade and 
thermal softening caused by increased temperature gradients. This 
blade, as is the case for most turbines, indicates that the thermal soften- 
ing effect outweighs the stiffening effect which results in lower frequencies 
for higher rotational speeds. 

The mode shapes from the analysis were normalized to give a unit gen- 
eralized mass for each mode. Mode shapes orthogonalized in this man- 
ner are referred to as normal modes. Table II shows the calculated 


natural frequencies for the three rotor speeds. 
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Plots showing the deformed blade mode shapes for these natural fre- 
quencies at 28,000 RPM are shown on figures 14 and 15. These modes 
can be classified using common plate theory modes as modes (1) first 
bending, (2) first edgewise, (3) first torsion and (4) second bending. The 
occurrence of the edgewise mode as the second weakest mode is unusual 
for turbomachinery. 

Plots of the mode shapes for strip no. 1, the strip closest to the blade tip 
are also shown on figure 16 as cross-sections through the blade airfoil 
section. These plots indicate that the first three modes are essentially 
rigid bending modes, with no appreciable chordwise bending of the 
airfoil. The fourth mode does display considerable bending of the airfoil 
section. These mode shapes displayed for this strip are the airfoil dis- 
placement descriptions which are enforced for each strip during the un- 
steady aerodynamic analysis. 

3.5 Aerodynamic Results 


The aerodynamic analysis for this blade was conducted for the 109% 
Rated Power Level (RPL) of the turbopump. This power level corresponds 
to a rotor speed of approximately 28,000 RPM. The SSME operates at this 
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power level for approximately 90% of the flight time, so the aeroelastic 
analysis is only presented for this rotor speed. 

The velocity triangles along the span of the blade were constructed based 
on information supplied by Rocketdyne. The turbine operates within a 
fully subsonic flow regime where the velocities represent a blade tip Mach 
number of 0.24 with a relative inlet Mach number of around 0.30. The 
working fluid is hydrogen-rich superheated steam, and it was assumed to 
obey the perfect gas law having a ratio of specific heats of y = 1.366. Ve- 
locity triangles for strip no. 1 showing the inlet and exit flows is included 
on figure 17. The superheated steam mixture has a sonic velocity at the 
turbine inlet of approximately a =5600 ft/sec and the inlet total pressure 
is P, = 4600 psi. 

The aeroelastic model for this turbine used the six strips indicated on 
figure 11, with the cascade properties used for each strip tabulated in ta- 
ble I. The steady potential flowfield was calculated for the strips sepa- 
rately and the computational grids and solutions stored in database files 
for use during the unsteady flow analysis. The computational grids used 
for the steady flow analysis consisted of a 78 x 25 mesh for the blade-to- 
blade global H-grid. A 70 x 1 1 C-grid was used for the local analysis grid, 


as described in section 2.11. 
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3.51 Steady Aerodynamic Results 

Contour plots showing the steady Mach number variation within the flow 
field for strips no. 1 and 3 are included as figures 18 and 19. The plots 
show the large expansion of the flow and acceleration around the suction 
surface of the blade, near the LE. The two plots are presented with the 
same contour level ranges, and very little variation in the steady flow field 
is seen between these two strips. 

Figure 20 shows the Mach number distributions along the surface of the 
airfoil. The figure includes the results for all six strips along the span of 
the airfoil. The large steady aerodynamic loading is evident for this tur- 
bine with the largest flow expansion near the leading edge of the suction 
surface. A similar figure indicating the airfoil steady surface pressure 
normalized by the inlet total pressure P r , is on figure 21. The surface 
pressures imply that the largest steady aerodynamic loading along the 
airfoil is in the outer 40 percent region of the airfoil near the tip. 

Inspection of the steady aerodynamic results in figures 20 and 21 shows 
that this turbine provides very little expansion and low increases in Mach 
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number from upstream to downstream, especially in comparison to air- 
craft gas turbines. The reason for this is that the dynamic pressure en- 
tering the blade is very high, approximately 480 psi at the tip. In addition, 
the inlet total pressure is close to 4600 psi, so that the static pressure 
nondimensionalized by total pressure can be misleading. The pressure 
difference across the airfoil can be up to 300 psi near the blade tip. These 
high pressures are a result of the superheated steam mixture having a 
high density and the high kinematic fluid velocities. 


3.52 Unsteady Aerodynamic Results 

The unsteady potential flows were calculated using the airfoil mode 
shapes from the finite element normal modes applied along each strip of 
the blade. At each strip, the flow reduced frequency parameter k was 
calculated based on the inlet relative flow velocity, strip airfoil chord 
length c and the assumed flutter frequency co f . The interpolated airfoil 
mode shape was used to prescribe the direction of airfoil oscillation as 
discussed in section 3.1. The unsteady aerodynamic program was used 
to determine the unsteady potential solution and unsteady harmonic 
pressures caused by the prescribed aeroelastic conditions. 
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The harmonic unsteady pressures resulting from blade motion do not re- 
spond in-phase with the blade oscillation. The phase lag or lead of the 
pressure is an important parameter which governs the stability of the 
blade. Unsteady pressure is represented as a complex quantity, with the 
real part representing the pressure which is in-phase and the imaginary 
part representing the pressure which is 90 degrees out-of-phase with the 
blade motion. 

The unsteady pressure distribution (real part) due to motion in the second 
mode, the edgewise mode, is shown along the blade span on figure 22. 
The ordinate is plotted as the unsteady pressure coefficient C p multiplied 
by the square root of the generalized mass for the mode. The in-phase 
pressure response varies strongly along the airfoil span, with the largest 
pressure gradients within about 30 percent of the leading edge. Note that 
the pressures are somewhat irregular, especially near the leading edge. 
This might be due to the discrete description of the airfoil shape. This 
apparent dependence upon airfoil description is a common problem 
among both steady potential and unsteady potential flow CFD analysis. 

The pressure distribution (imaginary part) for the same edgewise mode 
of motion is included in figure 23. The variation of the out-of-phase pres- 
sures, from the airfoil hub to the tip, obeys the same characteristics as the 
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real pressures, although these pressures are much smoother. A large 
change in the shape of the pressure distribution is seen between the last 
two strips along the blade span. 

For a single mode of motion, the stability of the unsteady flow can be de 
termined by inspecting the ratio of the real to the imaginary pressures at 
a particular location on the airfoil. If the pressures lag the blade motion, 
the motion is stable, if the pressures lead the motion it could become un 
stable. For systems which exhibit coupling between modes, this simple 
rule can no longer be used solely for determination of stability. An aero- 
dynamic work parameter (Verdon [38]) is used to determine the unsteady 
aerodynamic energy transfer for arbitrary modes of airfoil motion. The 
aerodynamic work per cycle represents a measure of the aerodynamic 
work provided by the fluid during one complete cycle of airfoil motion. A 
condition when the aerodynamic work per cycle is negative implies that 
the airfoil is performing work on the fluid. Conversely, a positive 
work/cycle signifies that the fluid is doing work on the airfoil, which can 
lead to an instability of the airfoil oscillation. 

Figure 24 shows the local work/cycle along the airfoil surface for the 
edgewise vibrational mode of the HPOTP. The work/cycle is unstable 
along almost the full chord length of the blade for all strips except the one 
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located at 58 percent of the airfoil span. These results would imply that 
the HPOTP turbine may be unstable in this edgewise mode of motion due 
to the motion-dependent response of the blade. 

3.6 Aeroelastic Results 


Calculations of the damped eigenvalues for the HPOTP first stage turbine 
were performed by solving the aeroelastic eigenvalue problem of 
equation (32). Four normal modes were retained within the modal analy- 
sis and the flutter analysis was performed when the assumed flutter fre- 
quency corresponded to the first three in-vacuum natural frequencies. 
The rotor was assumed to be tuned, and half of the total interblade phase 
angle modes were analyzed to reduce the computational time. The cal- 
culations were performed on the NASA LeRC Cray X-MP and Y-MP com- 
puter systems. The unsteady aerodynamic calculation, for a single strip 
and one interblade phase angle required approximately 12 CPU sec. The 
flutter analysis, for all strips at one interblade phase angle and one as- 
sumed frequency required 210 CPU secs. 

The first stability analyses were performed when the effect of mechanical 
damping due to the friction dampers was neglected. Flutter was assumed 
to occur at frequencies close to each of the first three natural frequencies. 
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The eigenvalue problem of equation (32) was solved assuming four basis 
normal modes, whereas the system equations were of size 4x4. Each of 
the interblade phase angle modes was solved individually, so that the 4 
x 4 system was solved 39 times for vibration at an assumed frequency. 
The solution was accomplished by assuming a flutter frequency co f , and 
the modal unsteady aerodynamic matrix was calculated followed by the 
eigensolution for the 4 eigenvalues. The aerodynamic damping for this 
turbine blade was extremely low, such that the damped natural frequency 
was nearly equal to the assumed frequency, therefore only one iteration 
was required for the flutter search. 

The calculations showed that the HPOTP blade was aeroelastically stable 
for all interblade phase angles when the vibration occurred at frequencies 
close to the first and third natural frequencies. An inspection of the cal- 
culated modal eigenvectors showed that the aeroelastic modes were al- 
most purely single-mode vibrations. The coupling due to the aerodynamic 
matrix [A] was not strong enough to cause any appreciable coupling in 
the modal space eigenvectors. 

The aeroelastic eigenvalues /i are plotted as a function of the interblade 
phase angle a in figure 25. This root locus is for the case when the vi- 
bration was at the first natural frequency, o),= 4748 Hz and the numbers 
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on this plot identify the different interblade phase angle modes for the 
rotor. Note that the eigenvalues are normalized by the assumed flutter 
frequency. Therefore the abscissa for this plot represents the system 
damping ratio and the ordinate represents the damped frequency ratio. 
The aerodynamic damping ratio for vibration at this mode varies from 
C = 0.16 percent of critical damping for a = 350.77° up to C = 0.40 percent 
for a = 161.54°. The results for vibration at the third normal mode re- 
sulted in lower damping ratios, although the rotor was stable for all phase 
angles (figure 26). 

Calculations when the vibration was assumed to occur at the second 
normal mode, the edgewise mode, showed that the rotor was unstable for 
almost all the interblade phase angles considered. The root locus plot for 
this edgewise motion of the blade is included in figure 27. These results 
indicate that it is possible for this blade to undergo flutter in an edgewise 
mode of vibration when mechanical damping from the friction dampers is 
neglected. The only stable modes are for interblade phase angles in the 
range from 332.31 deg to 346.15 deg. This instability is surprising be- 
cause the occurrence of flutter in an edgewise mode of vibration is rather 
uncommon. A combination of the high blade camber, large steady aero- 
dynamic loading and flow expansion, and the edgewise mode of vibration 
appear to lead to this form of instability. This instability would not be 
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identified by the methods cited in the literature review because of the use 
of the simplified aerodynamic and blade dynamic models. 


The potential for the HPOTP blades to experience flutter in this mode ap- 
pears to explain the occurrence of blade shank cracking noticed during 
the SSME development program. In the early 1980's, blade-to-blade fric- 
tion dampers were installed within the HPOTP rotor to attempt to damp 
out the unidentifiable vibrations. Continual high-cycle fatigue (HCF) 
cracking determined that a more efficient friction damper design was re- 
quired. A new two-piece damper was incorporated in the HPOTP design 
which reduced the cracking problem significantly. 

The mechanical damping provided by a blade-to-blade friction damper 
can be anywhere between 0.5 to 2.5 percent of critical damping (refs [43] 
and [44]). Measurements of the damper performance at Rocketdyne in a 
rotating dynamic test rig showed that the two-piece damper delivers al- 
most 1 percent of critical damping when the blades are stimulated in the 
edgewise vibrational mode. 

The stability calculations for the edgewise mode of vibration were contin- 
ued when the effect of the friction dampers was included. This entailed 
adding modal damping of { = 0.01 for the edgewise mode of the blade. 
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The modal damping ratios for the three other modes were assumed to be 
zero. The eigensolution was performed for the 39 interblade phase angle 
modes and the resulting root locus is shown in figure 28. It is observed 
that including the effect of modal damping essentially translates the root 
locus to the left (more stable) and slightly lowers the root locus (lower 
frequency). These results indicate that the HPOTP rotor is fully stable 
when the effect of the turbine blade friction dampers is included. The in- 
creased positive mechanical damping completely overwhelms the unsta- 
ble negative aerodynamic damping, resulting in the rotor becoming stable 
for all interblade phase angle modes. 

This result agrees with the history of the SSME HPOTP rotor, where early 
in the development program, blade cracking could have been due to ei- 
ther (a) flutter or (b) large forced response due to very low aerodynamic 
damping. When additional mechanical damping was introduced by using 
a more effective friction damper, the blade cracking problem was re- 
duced, a result in agreement with the analytical results presented here. 



CHAPTER 4 - CONCLUSION AND RECOMMENDATIONS 


4.1 Conclusion 

A modal aeroelastic analysis has been developed which has application 
to turbomachinery of general shape and over a range of flow regimes. 
The quasi-three-dimensional model combines the three-dimensional finite 
element modal results for the blade with stacked two-dimensional 
axisymmetric streamsurfaces along the span of the airfoil section of the 
blade. The steady flowfield along each strip was evaluated by solving for 
the nonlinear full potential flow using a finite-area solution method. 

The unsteady flowfield due to airfoil oscillation was determined using a 
linearized unsteady potential flow model. The unsteady pressures were 
integrated along the airfoil surface for each strip to result in modal un- 
steady aerodynamic forces. A modal expansion of the unsteady forces 
was used to determine the aerodynamic matrix which couples the modes 
in modal space. An iterative complex eigenvalue problem was expressed 
in state-vector form to include the effect of modal damping. 
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This modal aeroelastic approach was applied to the flutter analysis of a 

high-energy turbine blade from a rocket engine oxidizer turbopump. 

Specific results of this analysis were: 

• The natural frequencies for the turbine blade calculated using 
MSC/NASTRAN agreed very well with measurements reported by 
Rocketdyne. The modes occurred in the following order; (1) first 
bending, (2) first edgewise, (3) first torsion and (4) second bending. 

• Steady-state aerodynamic results indicated that this turbine is sub- 
jected to high steady gas loading, with pressure differences across 
the airfoil of up to 300 psi. The axial pressure drop through the 
blade-row from upstream to downstream was much smaller. 

• Unsteady aerodynamic results showed that the vibration of the blade 
at a single uncoupled mode, the first edgewise mode, would be un- 
stable because the fluid supplies energy to the oscillation which re- 
sults in a positive work of the fluid on the blade. 

• Modal flutter calculations determined that the tuned rotor was stable 
for vibration of the rotor at either the first or the third natural fre- 
quency, for all of the tuned interblade phase angle modes of the rotor. 
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• Results from the flutter computation for vibration at the second natural 
frequency indicated that flutter would occur for almost all of the 
interblade phase angle modes of the rotor. This result agrees with the 
high-cycle fatigue cracking problem encountered during SSME devel- 
opment in the early 1980's. 


• The addition of mechanical damping to simulate the blade friction 
dampers resulted in a stabilization of the self-excited vibration at the 
second natural frequency. 

4.2 Recommendations 


Further expansion of this model could include emphasis on the develop- 
ment of advanced unsteady aerodynamic models. In particular, models 
based on linearizations of rotational flows, i.e. linearized Euler solution, 
with emphasis on three-dimensional geometries is planned. The addition 
of a three-dimensional aerodynamic model will simplify the logistical 
problems associated with the aeroelastic strips and the modal aerodyna- 
mic force computation. The advantage of a linearization of the unsteady 
Euler equations is that unsteady flows containing shocks and rotational 
flows can be analyzed where the potential equation is no longer valid. In 
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addition, the computational cost of the linearized Euler solution is not 
prohibitive, although it is significantly larger than for the linearized po- 
tential solution. An example of an unsteady linearized Euler solution for 
two-dimensional turbomachinery flows was presented by Hall and 
Crawley [45], 

An effort is currently underway to analyze rotors which have small differ- 
ences in natural frequencies of the blades within the rotor, the so-called 
mistuned aeroelastic model. The mistuned rotor model is a more practi- 
cal representation of real-world turbomachinery since it takes into account 
the statistical differences in properties around the rotor. The occurrence 
of localized modes, where only a few blades respond with large ampli- 
tudes, can be captured using mistuning analyses. Blade failures typically 
occur with only a few blades cracking or fracturing, where these blades 
are referred to as "rogue" blade failures. A mistuned rotor analysis can 
provide better prediction of which blades may fail on a rotor due to flutter 
or forced response problems. 

Finally, a more accurate representation of the friction-damper model, in- 
corporating the micro-slip and macro-slip analyses could prove useful. 
The present method of accounting for friction-damping as equivalent 
modal damping is only a global representation of the effect of friction 
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dampers. Recent work has focused on applying nonlinear friction models 
using either lumped parameter or finite element and component mode 
models of blade elements. Some examples of the lumped parameter 
blade models using nonlinear friction damper elements are cited in Ref- 
erences [44] and [46]. 
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APPENDIX A 


DERIVATION OF LINEARIZED UNSTEADY POTENTIAL EQUATION 

The governing field equations for the flow of an inviscid, compressible 
fluid are described by the equations of mass, momentum, and energy 
conservation. These equations form the mathematical basis for de- 
scription of the Kinematic and thermodynamic modeling of fluid dynamics. 
The emphasis of the current work is to study the two dimensional, un- 
steady flow of an inviscid, compressible fluid. The equation set describing 
this flow is obtained by reducing the governing equations obtained from 
first principles. 

The governing equations for the unsteady, compressible, inviscid flow of 
a fluid are equations (A.1) to (A.3). 

— + V-(pV) = 0 (4.1) 

dt 

+ V VV = - \ Vp (4-2) 

dt P 

UL + v-Vh = + V-Vp (4.3) 

dt at 
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These equations are in conservation form and are general for a homoge- 
neous fluid which is flowing without viscous dissipation, no heat addition, 
and no body forces. 

The assumption that the fluid flows reversibly, and is thus isentropic, al- 
lows use of the perfect gas relationship which relates the local pressure 
to the density change as 


-^- = Constant (4.4) 

P 

The velocity potential is introduced by prescribing that the fluid flow is 
irrotational. This irrotationality condition is represented as 

Vxl/ = 0 (4.5) 

which permits defining the scalar velocity potential as 

V = (4.6) 

The equations (A.1-A.6) defined above represent the basis framework re- 
quired for the development of the steady-state and the unsteady potential 
equations utilized within this work. 
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UNSTEADY BERNOULLI-KELVIN THEOREM 

A relationship which defines the dependence of the fluid properties with 
the kinematic state is now required. This formulation is obtained from the 
momentum equation of equation (A.2). 

_|_ + vw = -i-Vp 

The following vector identity is useful for the reduction of this equation. 

(V-V)V =y V(V-V)- V x (Vx V) (A.7) 

Introduction of the irrotationality condition equation (A.5) and the velocity 
potential definition equation (A.6) reduces the above vector relationship 
to 

(\/ ; V)\7 = yV(Vtf>V</>) (A.8) 

Substitution of the above equation into the momentum equation (A.2), re- 
sults in 

(A. 9) 


■jf (v<« + y V(V0-V0) = -j- v P 
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V0 f + V 


(v<W) 


+ — Vp = 0 


(/\. 10 ) 


Now, recalling from the calculus, since 


fVp = V 



(A. 11) 


substituting equation (A.11) results in the following integral equation 


, ( v<f>-V4> ) 

V4> t + v( ■ — ] + v 


(K)- 


(A. 12) 




0 


and integrating gives 


(V0-V0) 


(Pt + 2-I_Ll + jv^- = G(f) 


(A. 13) 


This equation is known as the Bernoulli-Kelvin equation and it applies 
along a streamline in the fluid. 
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UNSTEADY FULL POTENTIAL EQUATION 

The unsteady potential equation is derived by expansion and simplifi- 
cation of the mass conservation equation indicated in equation (A.1). This 
equation is reproduced below for completeness. 

— + V(pV) = 0 ( /A - 1 ) 

dt 

The divergence term of this equation may be expanded and the equation 
rearranged to result in 


+ v v + J-\/Vp = o ( yA - 14 ) 

p dt p 


The first term of this equation may be written as the product of two partial 


derivatives using the chain rule 


j_ 1 d P d P (4.15) 

P dt P dp dt 


and assuming the fluid behaves as a perfect gas, under isentropic condi- 
tions, the sonic velocity becomes 


a 


dp_ 

dp 


s = constant 


(A. 16) 
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such that equation (A. 15) now becomes 


1 _ = 1 dp 

P dt ~ pa 2 dt 


(A. 17) 


cp 

A relationship for is determined through the use of the Bernoulli- 
Kelvin Theorem of equation (A. 13). Differentiating equation (A. 13) with 
respect to time results in 


$tt + 


d_ 

dt 


2 ) dt J p 


(A. 18) 


where the derivative of the time constant G(t) has arbitrarily been set to 
zero. The last term of equation (A. 18) can be written as 


d_ f y _P 1 dp 

dt J P P dt 


(A. 19) 


and so rearranging equation (A. 18) and substituting equation (A. 19) gives 


dp_ 

dt 



v<W ^ 


(A. 20) 


and substituting equation (A. 20) into equation (A. 17) results in the follow- 
ing relation for the first term in equation (A. 14) 
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1 d£ 

P dt 


- -V + vw,) 

a 


(/A .21) 


The second term from equation (A.14) is simply the Laplacian of the po- 
tential 

V-\/=V 2 (/> (A .22) 


The last term of equation (A.14) can be written as 


V-Vp Vtf>Vp 


pa 

V0-V 


K 


(A. 23) 


solving equation (A.12) for V (/^T") and sustituting into equation 


and 
(A. 22) gives 


V-Vp V0-V / , 

-s-= — + — r~ 


{A. 24) 


At this point, each of the individual terms from the continuity equation of 
equation (A.14) have been expanded. These terms of equation (A. 21), 
(A. 22), (A. 24) are substituted back into equation (A.14). 


4 U + 2 V<W, + -i- V*-V(V0-V0) = a 2 V 2 ,#. 


(4.25) 


This equation is known as the unsteady full potential equation. It contains 
only terms of the scalar velocity potential 0 and the local fluid sonic ve- 
locity a. 


An expression for the sonic velocity can be obtained from the Bernoulli- 
Kelvin equation developed previously. Along a streamline, from a point 
at far upstream infinity where the flow is fully steady, the Bernoulli 
equation (A. 13) becomes 


a oo + (y- *) 


1/ 


4>t - -5- W-v</>) 


(A .26) 


So the above equations (A.25) and (A. 26) completely describe the un- 
steady potential flow for a compressible) inviscid fluid. 


LINEARIZED POTENTIAL EQUATIONS 


The time varying potential is assumed to be a series expansion in time of 
harmonic terms as shown below 

0(x,y,f) = d>(x,y) + (f)(x ,y)e lcot + 0(£ 2 ) (4-27) 

where c is of the order of the small perturbations r, and </>(x,y). This as- 
sumption results in a zeroth-order term which represents the steady-state 
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potential $(x,y), and a first-order harmonic unsteady potential 
</)(x,y,f)e ,cu '. The terms of higher than first order in e, i.e. higher than e' Mt , 
are truncated from the expansion. This expansion can now be substituted 
into the unsteady full potential equation (A.25). 

The following expressions will be useful for this development. They are 
obtained by applying the expansion of equation (A.27). 


4> tt = - co 2 <t> 
y<f) = VO + Vtf> 

V</>, = /wVj ( yA - 28 ) 

v 2 </> = v 2 a> + v 2 4> 

V0-V</> = (VO) 2 + (vj) 2 + 2 vo-vj 

The above relations can be used to expand the unsteady full potential 
equation to result in an equation which contains terms of both 0(0) and 
order O(e). These terms are grouped according to their order into two 

equations. 

The equation consisting of terms of order 0(0) is known as the steady full 
potential equation as shown below. 


-1 VO V(V<D VO) = a 2 V 2 0 


(A .29) 
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The same equation can be obtained by dropping the time-dependent 
terms from the unsteady full potential equation of equation (A.25). 

An equivalent equation for the terms of order 0(f) is shown below 

- c o 2 4> + 2/coVO-vJ + -y V0-V(V<D) 2 + V<t>-V(V$vJ) 

1 {A .30) 

4- vJ-V(VO-V<£) = a 2 V 2 J 

This equation represents the linearized unsteady potential equation, 
where the potential $(x,y) is assumed to be harmonic in time. 

The substantial derivative operator for the potential having harmonic time 
dependence becomes 


D<{) 

~Dt 


— icocj) + VOV</> 



—co 2 ({) + 2/'coV<t>V</> + V<D-V(V<I>-V</>) 
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The substitution of the above substantial derivative operators into the un 
steady potential equation (A.30) simplifies to 


2 ^ ~ 

D Z - + -! vJ-V(Vd)) 2 + vj>-V(V<l>V</>) = a 2 V 2 </> (A.31) 

Dt 2 2 


This equation represents the linearized unsteady potential flow where the 
unsteady potential 4> is assumed to be harmonic in time and of small am- 
plitude. 
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Figure 2. Cascade and airfoil geometry and nomenclature 
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(a) 



Figure 4. Steady flow computational meshes, (a) global mesh, 
(b) local mesh 
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Undeformed Airfoil 



Figure 6. Oscillating airfoil displacement vector r definition 
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Figure 7. FREPS program aeroelastic stability logic flowchart 
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Figure 8. Modal aerodynamic matrix [/A] logic flowchart 
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Figure 9. SSME high pressure oxygen turbopump (HPOTP) cross-section 
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Figure 10. SSME HPOTP turbine blade finite element model 
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Figure 11. Turbine blade aeroelastic strip definitions 
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Figure 12. Turbine blade airfoil description at 75% span 
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Figure 13. Calculated natural frequencies versus rotational speed 





Figure 14. Calculated deformed mode shapes for modes 1 and 2 



Figure 15. Calculated deformed mode shapes for modes 3 and 4 
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Figure 16. Airfoil mode shape cross-sections at strip no. 1 
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Figure 17. Turbine blade fluid velocity diagram at strip no. 1 
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Figure 18. Steady flow Mach number contours for strip no. 1 
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Figure 19. Steady flow Mach number contours for strip no. 
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Figure 23. Turbine blade unsteady surface pressure distributions 
(imaginary part) due to motion in second mode (edgewise) 




Figure 25. HPOTP turbine rotor tuned root locus for vibration 
in first mode (bending) 
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Figure 27. HPOTP turbine rotor tuned root locus for vibration in 

second mode (edgewise), mechanical damping neglected 
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Figure 28. HPOTP turbine rotor tuned root locus for vibration in 

second mode (edgewise), 1.0% modal damping ratio included 
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